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Monte  Carlo  simulations  are  used  to  model  the  self-organizing  behavior  of  the  biomineralizing 
peptide  KSL  (KKVVFKVKFK)  in  the  presence  of  phosphate.  Originally  identified  as  an 
antimicrobial  peptide,  KSL  also  directs  the  formation  of  biosilica  through  a  hypothetical 
supramolecular  template  that  requires  phosphate  for  assembly.  Specificity  of  each  residue  and  the 
interactions  between  the  peptide  and  phosphate  are  considered  in  a  coarse-grained  model. 

Both  local  and  global  physical  quantities  are  calculated  as  the  constituents  execute  their  stochastic 
motion  in  the  presence  and  absence  of  phosphate.  Ordered  peptide  aggregates  develop  after 
simulations  reach  thermodynamic  equilibrium,  wherein  phosphates  form  bridging  ligands  with 
lysines  and  are  found  interdigitated  between  peptide  molecules.  Results  demonstrate  that 
interactions  between  the  lysines  and  phosphate  drive  self-organization  into  lower  energy 
conformations  of  interconnected  peptide  scaffolds  that  resemble  the  supramolecular  structures  of 
polypeptide-  and  polyamine-mediated  silica  condensation  systems.  Furthermore,  the  specific 
phosphate-peptide  organization  appears  to  mimic  the  zwitterionic  structure  of  native  silaffins 
(scaffold  proteins  of  diatom  shells),  suggesting  a  similar  template  organization  for  silica  deposition 
between  the  in  vitro  KSL  and  silaffin  systems. 


1.  Introduction 

The  mechanism  of  silica  formation  has  been  long  studied  in 
the  materials  sciences,  where  the  polycondensation  of  silicic 
acid  can  be  completed  using  acid-base  chemistry  and  directed 
into  a  variety  of  silica  structures  using  organic  molecule 
mediators.1  Similar  condensation  chemistry  occurs  in  natural 
systems  at  the  nanoscale  and  is  dominated  by  marine  micro¬ 
organisms  such  as  diatoms,  radiolarian,  and  synurophytes.2,3 
The  microbes  perfected  silica  synthesis  millions  of  years  before 
us  and  to  a  level  of  complexity  that  is  superior  to  our  ability  to 
create  similar  structures  in  the  laboratory.  For  the  materials 
chemist,  the  goal  is  to  mimic  the  architecture  of  these  biological 
systems  for  useful  applications  in  nanotechnology. 
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f  Electronic  supplementary  information  (ESI)  available:  Larger  snap¬ 
shot  images  at  t  =  105  and  106  for  simulations  containing  peptide 
concentrations  0.01,  0.02,  0.03,  0.05,  and  0.07.  See  DOI:  10.1039/ 
c0cp01364a 


In  the  past  two  decades,  research  unlocked  many  mysteries 
of  diatom  shell  (frustle)  formation.  Critical  discoveries  re¬ 
vealed  that  the  polyionic  proteins  (silaffins)  and  long-chain 
polyamines  (LCPAs)  synthesized  by  select  diatom  species  are 
used  to  precipitate  silica  and  direct  its  deposition  into  exquisite 
architectures.4^6  When  extracted  from  the  frustules  of  diatoms, 
silaffins  promote  in  vitro  silica  precipitation  with  silicic  acid  in 
aqueous  solutions  at  a  wide  pH  range.7  Silica  forms  rapidly  at 
room  temperature  in  the  absence  of  other  ionic  catalysts  or 
cofactors  (buffers,  salts,  and  other  ionic  species)  that  are 
normally  required  with  in  vitro  biosilica  synthesis.  The  silaffins 
are  uniquely  functionalized;  serine  residues  are  phosphorylated 
and  select  lysines  are  posttranslationally  modified  with 
branched  poly  amine  chains.  This  zwitterionic  structure  of 
the  protein  facilitates  self-assembly  into  macromolecular 
structures,  which  mediates  the  polycondensation  of  silicic  acid 
to  form  silica  spheres  in  vitro  and  intricate  diatom  frustules 
in  vivo. 

Since  the  discovery  of  silaffins  and  LCPAs,  several  other 
unrelated  proteins,  synthetic  peptides,  and  polyamines  have 
also  been  shown  to  catalyze  silica  biomineralization  reactions.8,9 
Synthetic  and  other  natural,  unmodified  peptides  that  induce 
silica  condensation  share  similar  electrostatic  properties  to 
natural  silaffins,  but  most  require  the  addition  of  multivalent 
anions  ( e.g .,  phosphate)  and  are  restricted  to  circumneutral 
pH  to  facilitate  biomineralization  reactions.10  12  The  presence 
of  phosphate  and  a  combination  of  charged  and  uncharged 
amines  in  neutral  to  alkaline  conditions  mimic  the  properties 
of  the  poly-ionic  modifications  on  native  silaffins.13,14  The 
zwitterionic  properties  facilitate  intermolecular  interactions 
and  promote  self-assembly  into  supramolecular  aggregates, 
which  has  been  shown  to  be  a  prerequisite  for  silica  formation. 
Liquid-phase  31P  NMR  spectroscopy  of  the  native  silaffin 
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1A  peptide  (from  the  frustle  of  the  diatom  Cylindrotheca 
fusiformis)  and  dynamic  light  scattering  analysis  on  polyamines 
demonstrated  that  these  silica-forming  molecules  are  not 
monomeric  in  solution  and  form  aggregates.6,10,11,13  Other 
studies  suggest  that  self-assembly  is  required  for  biosilica 
formation,  as  it  provides  a  structure- specific  template  for  rapid 
silicic  acid  poly-condensation  or  forms  microemulsions  that 
act  as  flocculating  agents  to  promote  silica  aggregation  and 
precipitation  via  a  phase  separation  mechanism.12,15,16  In  a 
seminal  study  by  Kroger  et  al .,  dephosphorylated  silaffins  did 
not  self-assemble  and  thus,  had  lost  the  ability  to  biomineralize 
silica.7  The  addition  of  phosphate  ions  restores  silica  formation 
by  allowing  the  positively-charged  peptides  to  self-assemble 
via  intermolecular  phosphate  anion  bridges.  Likewise,  the 
unphosphorylated  peptides  that  mimic  silaffins  require  the 
addition  of  phosphate  to  catalyze  silica  formation.  Under¬ 
standing  the  mechanism  of  self-assembly  prior  to  silica 
condensation  is  the  first  step  in  synthesizing  complex  silica 
structures  at  the  molecular  level. 

The  antimicrobial  peptide  KSL  (KKVVFKVKFK)  is  a 
pertinent  example  of  a  peptide  that  catalyzes  fortuitous 
biomineralizations  reactions.17  KSL  will  promote  silica 
precipitation  and  form  silica  nanoparticles  when  added  to 
solutions  containing  phosphate  buffer  (pH  =  8)  and  hydrolyzed 
tetramethyl  orthosilicate  (TMOS).  The  sequence  of  KSL  is 
rich  in  lysines  that  are  interspersed  with  the  hydrophobic 
amino  acids,  valine  and  phenylalanine.  In  addition  to  its 
biomineralization  properties,  KSL  is  a  potent  antimicrobial 
agent.  The  sequence  endows  the  peptide  with  specific  ionic  and 
hydrophobic  properties  that  imparts  biocidal  activity  toward  a 
wide  range  of  microorganisms.18  Cationicity  and  amphiphilicity 
are  inherent  to  many  antimicrobial  peptides  (AMPs)  and 
several  models  suggest  that  AMPs  will  self-assemble  into 
supramolecular  structures  within  microbial  cell  membranes 
and  disrupt  their  integrity  and  function.18,19  While  the  anti¬ 
microbial  activity  and  biomineralization  properties  of  KSL  are 
likely  coincidental,  we  propose  that  the  self-organization  due 
to  cationic  and  amphiphilic  properties  of  KSL  also  mediates 
rapid  silica  condensation. 

Monte  Carlo-led  lattice  simulations  have  been  extensively 
used  to  describe  the  aggregation  properties  of  proteins  and 
small  peptides  in  amyloid  and  prion  diseases,20,21  interactions 
between  AMPs  and  cellular  membranes,22,23  and  bio-inorganic 
interactions  that  occur  in  peptide-mediated  biomineralization 
reactions.24^28  Generic  aspects  of  oligomerization  and  patterned 
assembly  can  be  gleaned  from  simulations  that  provide  simplified 
protein  models  based  on  relevant  properties  of  protein-protein 
or  protein-solvent  interactions  and  computational  algorithms 
that  rely  on  repeated  random  sampling  of  stochastic  motion. 
In  a  previous  study,  Lenoci  and  Camp  have  analyzed  the  self- 
assembly  of  silaffins  peptides  with  long-chain  polyamines  by  a 
coarse-grained  computer  simulation  model.24  Using  an 
attenuated  description  of  polyamine  tails  on  the  native  silaffin 
1A  protein,  the  self-organization  of  peptides  was  successfully 
captured  by  coarse-grained  interactions  among  different 
components.  The  Lenoci  and  Camp  model  varied  the  concen¬ 
tration  of  the  peptide  fragment  and  studied  its  aggregation 
leading  to  scaffolds  (clusters)  and  complementary  pore  spaces 
(cavities).  The  peptide  clusters  are  isolated  with  open  pores  at 


low  concentration  and  large  connected  clusters  with  confined 
cavities  at  higher  concentrations.  The  results  suggest  that 
self-assembled  structures  act  as  either  nucleation  points  or 
scaffolds  for  the  deposition  of  silica. 

In  our  study,  we  consider  the  first  computational  model  to 
understand  the  self-organization  of  a  biomineralizing  peptide 
in  the  presence  of  an  ionic  solution.  Unlike  the  study  by 
Lenoci  and  Camp,  we  retain  the  electrostatic  and  hydrophobic 
properties  of  each  amino  acid  and  incorporate  explicit  solvent 
to  represent  phosphate.  While  the  structural  details  of  each 
amino  acid  are  ignored,  the  model  incorporates  the  specific 
strength  of  every  possible  interaction  between  two  different 
amino  acids,  based  upon  a  database  of  known  interactions  (see 
below).  Because  of  the  efficiency  of  the  bond-fluctuation 
model29,30  with  ample  degrees  of  freedom  on  a  discrete  lattice, 
we  are  able  to  explore  a  relatively  larger  phase  space  by  the 
Monte  Carlo  simulations  to  study  the  self-organizing  structural 
evolution  with  a  range  of  interaction  strengths  between 
phosphate  and  lysine.  We  examine  a  number  of  local  and 
physical  quantities.  The  results  provide  insight  into  the 
self-organizing  behavior  of  short,  amphipathic  peptides  in 
biological  systems  and  in  nanomaterials  synthesis. 


2.  Model 


A  coarse-grained  description  for  KSL  is  used  to  examine  the 
self-organizing  structures  in  the  presence  and  absence  of 
phosphate  on  a  three  dimensional  discrete  (cubic)  lattice.26,31 
An  amino  acid  residue  is  represented  by  a  particle  or  node. 
Therefore,  the  peptide  is  composed  of  a  set  of  nodes  tethered 
together  in  a  flexible  chain  in  a  specific  sequence.  KSL  is  a 
chain  of  ten  amino  acids:  five  lysines  (Lys,  K),  three  valines 
(Val,  V),  and  two  phenylalanines  (Phe,  F),  in  the  specific 
sequence,  KKVVFKVKFK.  The  peptide  chain  with  the 
covalently  bonded  nodes  is  represented  by  the  bond-fluctuation 
model  on  a  lattice  with  a  cubic  grid.29  A  node  is  represented  by 
a  unit  cube  (occupying  its  eight  lattice  sites)  and  the  bond 
length  between  consecutive  nodes  can  fluctuate  between  2  and 

10  with  an  exception  of  ^8  in  unit  of  lattice  constant.  Such 
a  bond-fluctuation  description  is  known  to  capture  the 
computational  efficiency  while  incorporating  ample  degrees 
of  freedom  in  complex  polymer  systems,  peptides,  and  multi- 
component  nanocomposites.26,29,31,32  The  phosphate  constituent 
is  represented  by  a  particle  with  the  same  size  as  that  of  a 
residue  node. 

The  simulation  box  consists  of  peptides  and  phosphate 
constituents  with  concentration  CKsl  and  CPh,  respectively. 
Apart  from  excluded  volume  constraints,  each  residue  and 
phosphate  constituent  interact  with  a  pair  potential.  The 
interaction  strength  between  two  elements  at  sites  i  and  j 
separated  by  a  distance  ry  is  represented  by  the  Lennard-Jones 
(LJ)  potential: 


Uij=f 


rij<r  c 


where  rc  ( =  yj 8)  is  the  interaction  cut-off  range  and  a  =  1  in 
unit  of  lattice  constant.  The  strength  of  the  interaction  fay  is 
controlled  by  an  arbitrary  parameter  f  (f  >  0)  that  can  be 
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varied  with  positive  (repulsive)  and  negative  (attractive) 
parameter  Sy  characteristics  of  each  constituent  (see  below). 
The  distance  between  two  sites  is  discrete  and  measured  in 
lattice  constant  units.  The  depth  of  the  LJ  potential  is 
controlled  by  the  magnitude  of  the  pair  interaction  strength 
f&ij  in  Boltzmann  constant  units.  The  interaction  matrix  could 
be  relatively  large  for  large  peptides.  Since  only  three  different 
amino  acid  residues  (K,  V,  F)  comprise  KSL,  the  interaction 
matrix  is  relatively  small.  The  number  of  the  pair  interaction 
matrix  elements  is  further  reduced  due  to  symmetries 
(i e.g .,  Sy  =  £jt).  The  magnitude  of  each  interaction  element 
(sy)  is  based  on  the  insight  gained  from  all  atomistic 
description  and  general  characteristics.26  To  capture  the 
specificity,  every  possible  interaction  between  two  different 
amino  acids  is  considered.  The  interaction  strength  of  each 
residue-residue  contact  is  based  on  a  table  of  inter-residue 
contact  energies  derived  from  a  set  of  known  residue  inter¬ 
actions  in  the  Protein  Data  Bank.  In  this  study  by  Miyazawa 
and  Jernigan,33  inter-residue  contact  energies  for  all  pairs  of 
amino  acids  were  calculated  empirically  by  quantifying  and 
averaging  all  residue-residue  contacts  (Nrr  >  18  000)  identified 
in  a  set  of  globular  protein  crystal  structures.  Based  on  the 
assumption  that  the  native  conformations  of  proteins  reflect 
the  lowest  free  energy  of  residue-residue  interactions,  the 
average  of  all  inter-residue  contacts  found  in  the  population 
of  protein  structures  was  used  to  formulate  a  relative  table  of 
intrinsic  interaction  strengths  for  all  amino  acid  pairs. 

In  KSL,  most  residues  are  lysine  (K)  and  are  positively 
charged,  while  the  remaining  residues  are  phenylalanine  (F) 
and  valine  (V)  and  are  both  hydrophobic.  Based  on  the 
inter-residue  contact  energies  established  by  Miyazawa  and 
Jernigan,33  the  interactions  among  hydrophobic  groups  are 
attractive  and  their  interaction  strengths  are:  eVF  =  —0.22, 
£vv  =  —0.29,  and  eFF  =  —0.44.  Interactions  between  lysines 
are  repulsive  (sKk  =  +  0.25),  as  are  the  interactions  between 
lysine  and  the  hydrophobic  residues  (eVK  =  £Fk  =  +0.44). 
The  interactions  between  lysine  and  phosphate  are  more 
attractive  and  are  varied  in  simulations  (gKph  =  -0.5  to 
—0.8).  The  interaction  strength  between  lysine  and  phosphate 
is  varied  because  an  intrinsic  energy  of  this  pair  was  not 
determined  in  the  Miyazawa  and  Jernigan  study.33  Accordingly, 
we  choose  a  range  that  is  within  the  range  of  contact  energies 
calculated  between  oppositely-charged  amino  acids  (i.e., 
aspartic  and  glutamic  acid  vs.  lysine,  arginine,  and  histidine), 
which  ranged  from  approximately  —0.4  to  —1.0.  Phosphate 
buffer  at  physiological  conditions  is  mostly  either  H2P04_  or 
HPO4-2,  and  the  polyionic  nature  of  these  species  invokes 
transient  bonding  with  each  other,  so  phosphate-phosphate 
interactions  are  assigned  as  weakly  attractive  (ePhPh  =  -0.1). 
The  arbitrary  parameter  /  can  be  selected  to  accentuate  the 
differences  in  dependence  of  the  physical  quantities  on  the 
peptide  concentration;  we  choose  /  =  100  here.  The  reduced 
temperature  T  =  1  in  unit  of  the  Boltzmann  constant. 

Each  simulation  is  initiated  with  a  random  distribution  of 
peptide  chains  and  phosphate  particles  on  a  cubic  lattice  (L3). 
The  Metropolis  algorithm  is  used  to  stochastically  select  and 
move  each  peptide  node  and  phosphate  particle.  A  peptide 
node  or  phosphate  particle  at  a  site  i  and  one  of  its  neighboring 
site  j  are  selected  randomly.  If  the  neighboring  site  j  is  empty 


and  the  selected  constituent  at  site  i  is  a  phosphate  particle, 
then  an  attempt  is  made  to  move  the  phosphate  from  site  i  to 
site  j  with  the  Boltzmann  probability  e~AE/T,  where  A E  =  Ej  —  Et 
is  the  change  in  energy  as  a  result  of  the  move  and  T  is  the 
temperature  in  Boltzmann  constant  units.  If  the  selected 
constituent  is  a  peptide  node,  then  the  bond-fluctuation 
constraint  is  implemented  in  addition  to  Boltzmann  probability 
for  its  move  from  site  i  to  site  j.  A  time  step  (t)  is  defined  as  an 
attempt  to  move  each  residue  and  phosphate  particle  once 
(i.e.,  Monte  Carlo  step  or  MCS). 

Simulations  with  KSL  were  first  performed  in  the  absence  of 
phosphate.  A  series  of  simulations  were  then  carried  out  with 
varying  volume  fractions  of  peptide  chains  (CKsl  =  0.01-0.07) 
with  a  constant  phosphate  particle  volume  fraction  (CPh  =  0.1) 
for  a  range  of  interaction  strengths  between  lysine  and  phos¬ 
phate.  In  laboratory  experiments,  the  molar  ratio  between 
phosphate  and  KSL  is  much  higher  than  in  our  simulations 
(~  105).  In  order  to  study  peptide  aggregation  and  scaffolding, 
an  ample  number  of  KSL  chains  is  required  in  the  lattice,  in 
order  to  provide  a  sufficient  number  of  peptide  interactions  for 
meaningful  analysis.  At  experimental  molar  ratios,  the 
simulation  would  not  contain  an  adequate  number  of  peptide 
constituents.  To  overcome  the  deficiency,  more  peptide  chains 
can  be  included  in  a  larger  lattice  (e.g.,  1283  or  greater)  and 
simulation  times  extended  to  allow  the  peptides  to  interact  and 
reach  equilibrium.  Unfortunately,  the  required  increase  in 
lattice  space  and/or  time  steps  to  accommodate  experimental 
conditions  is  constrained  by  computing  hardware,  resources, 
and  time.  After  running  several  test  simulations  at  different 
lattice  sizes  (643  and  1283)  and  different  time  steps  (t  =  105  and 
106  MCS),  ultimately,  a  lattice  size  of  643  and  a  time  step 
length  of  t  =  106  were  selected  for  simulations  in  the  study. 
We  did  not  observe  significant  differences  in  peptide  dynamics 
between  simulations  using  different  lattice  sizes,  but  did  notice 
that  simulations  required  >  105  time  steps  to  reach  equilibrium 
(Fig.  SI  of  the  ESIf).  Results  are  an  average  of  10-50 
independent  simulations  for  each  sample  condition.  Under 
these  conditions,  the  consolidated  results  were  sufficient  to 
calculate  the  average  values  of  a  range  of  local  and  global 
physical  quantities  and  show  clear  trends  between  the  different 
simulation  conditions. 

3.  Results  and  discussion 

3.1  Self-assembly  of  peptide  in  the  absence  and  presence  of 
phosphate 

The  final  snapshots  (i.e.,  at  the  end  of  million  time  steps)  of 
simulations  containing  peptide  in  the  absence  and  presence  of 
phosphate  are  presented  in  Fig.  1.  In  the  absence  of 
phosphate,  a  minimal  amount  of  peptide  aggregation  is  observed, 
particularly  at  the  larger  CKsl  concentrations.  As  the  only 
attractive  forces  in  these  simulations,  it  is  likely  that  the 
hydrophobic  residues  interact  and  cause  the  peptides  to  relax 
into  small  aggregates.  The  peptide  chains  exhibit  mostly 
random  conformations,  although  a  few  are  shown  folded  upon 
themselves.  There  is  a  greater  probability  that  hydrophobic 
residues  on  the  same  peptide  will  interact  first,  resulting  in  the 
folded  conformation. 
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Fig.  1  Final  snapshots  (/  =  106  MCS)  for  simulations  of  increasing  peptide  concentration  (Cksl,  rows  top  to  bottom)  and  increasing 
lysine-phosphate  interaction  strength  (eKPh,  second  column  from  left  to  right).  First  column  on  left  is  final  snapshots  of  simulations  without 
phosphate.  All  snapshots  (including  CKsl  =  0-01  and  0.07)  are  provided  as  larger  figures  in  Fig.  S2  of  the  ESI.f 


In  simulations  containing  phosphate,  the  images  clearly 
show  that  phosphate  molecules  induce  peptide  to  organize 
into  clusters.  Large  aggregates  are  observed  in  snapshots  at  the 
lowest  lysine-phosphate  interaction  strength  (eKph  =  —0.5). 
Surprisingly,  increasing  the  eKph  appears  to  disrupt  the  local 
density  of  the  large  aggregates,  rather  than  causing  increased 
density  and  further  expansion  of  large  aggregates.  The  size  of 
the  aggregate  particles  decreases  and  the  network  of  peptides 
also  spreads  or  fractionates  into  smaller  complexes.  The 
competition  between  the  residue-residue  and  phosphate- 
lysine  interaction  leads  to  distinct  assemblies  of  aggregates  as 
a  function  of  interaction  strength.  At  low  gp-Ph,  hydrophobic 
interactions  play  a  role  in  aggregation  kinetics  and  yield  one 
type  of  aggregate  morphology.  At  high  eKPh,  it  appears  that 
phosphate-lysine  interactions  overcome  the  hydrophobic 
interaction  strength,  which  results  in  a  different  aggregate 
structure.  This  assessment  is  based  on  visual  interpretation 
of  the  simulation  snapshots;  in  the  following  sections,  we 
determine  the  quantitative  structural  correlations. 

3.2  Local  structural  correlations 

Correlation  between  local  structure  and  the  number  of 
interactions  between  different  constituents  is  determined  by 
calculating  the  average  number  of  hydrophobic  residues  (7VH) 
and  phosphates  (/VPh)  within  the  range  of  interaction  (r  =  ^8) 
of  each  residue.  In  simulations  containing  CKsl  =  0-03  and 
0.054  we  see  at  £KPh  =  —0.5,  the  hydrophobic  residues  of  KSL 
(V3,  V4,  F5,  V7,  and  F9)  contain  the  highest  number  of 

i  Unless  otherwise  shown,  results  from  simulations  containing 
Cksl  =  0.03  and  0.05  were  representative  of  simulations  at  all  CKsl- 


hydrophobic  interactions  (Fig.  2A).  As  eKph  increases,  AH  of 
the  hydrophobic  residues  decreases,  while  the  NH  around 
lysine  residues  (K1?  K2,  K6,  K8,  and  K10)  remains  relatively 
the  same.  The  opposite  trend  is  true  when  we  calculate  the 
average  number  of  phosphates  (NFh)  within  the  range  of 
interaction  around  each  residue  (Fig.  2B).  As  eKph  increases, 
the  number  of  phosphate  interacting  with  lysine  increases, 
while  the  number  of  phosphates  surrounding  the  hydrophobic 
residues  remains  relatively  constant.  This  trend  can  be  further 
clarified  by  plotting  the  average  number  of  hydrophobic 
residues  and  phosphates  within  the  range  of  interaction 
for  each  residue  as  a  function  of  eKph  (Fig.  3).  At  low 
^KPh?  hydrophobic  interactions  (labeled  (a)  in  Fig.  3)  are 
nearly  as  numerous  as  lysine-phosphate  interactions  (b). 
As  sKph  increases,  phosphate-lysine  interactions  dominate 
and  consequently,  NH  of  Val  and  Phe  decreases  (c).  Concurrently, 
the  number  of  phosphates  interacting  with  the  peptide 
increases  (d).  Here,  it  is  clear  that  the  aggregation  properties 
at  low  eKph  are  influenced  by  hydrophobic  interactions  and  at 
larger  f:KPh,  phosphate  percolates  into  the  clusters  and  re-organizes 
the  peptide  molecules  to  form  new  clusters  of  a  smaller  size 
and  lower  local  peptide  density. 

3.3  Distribution  of  constituents  within  peptide  aggregates 

Although  the  local  structural  correlations  discussed  above  are 
insightful,  it  does  not  provide  the  overall  structure  of  the 
whole  system.  To  gain  insight  into  the  spatial  distribution  of 
substituents  throughout  the  entire  lattice,  the  radial  distribution 
function  ( rdf)  can  be  used  to  correlate  the  global  structure  of 
the  system.  In  our  simulations,  the  rdf  measures  how  the 
density  of  residues  and  phosphates  changes  as  a  function  of 
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Fig.  2  The  average  number  of  hydrophobic  residues  (A)  and  phos¬ 
phates  (B)  around  each  residue.  Peptide  concentrations  in  simulations 
(Cksl)  are  0.03  (left)  and  0.05  (right)  and  the  lysine-phosphate 
interaction  strength  (eKPh)  for  each  bar  graph  is  shown  in  parentheses. 
Each  residue  is  represented  by  a  number  (shown  at  bottom). 


•  Nh  of  val  and  phe  ▲  NPh  of  val  and  phe 

■  Nh  of  lys  ▼  NPh  of  lys 


Fig.  3  Average  change  in  VH  and  NP h  for  hydrophobic  residues  and 
lysines  as  a  function  of  eKph-  Letter  designations  (a-d)  are  described  in 
the  text.  Simulation  conditions  are  the  same  as  described  in  Fig.  2  and 
error  bars  show  standard  error  of  the  mean. 

the  distance  from  each  substituent  (residue  or  phosphate).  The 
number  of  lattice  sites  either  occupied  by  a  residue  or  a 
phosphate  (Ah)  is  counted  between  two  concentric  spheres 


having  radii  rx  and  r2  (where  r2  >  rx)  measured  from  each 
substituent.  The  average  density  p(r )  of  this  region  is  Nx/N, 
where  N  is  the  total  number  of  lattice  sites  (empty  or  occupied) 
between  rx  and  r2.  We  can  evaluate  p(r)  around  each  residue 
within  a  range  of  r  lying  between  rx  and  r2  while  keeping  Ar 
sufficiently  small  and  estimate  the  average  (/.<?.,  rdf).  We  have 
used  r2  =  1-100,  which  allows  the  function  to  extend  further 
than  the  local  structural  correlation  and  gauge  the  distribution 
of  different  components  throughout  the  aggregates. 

Fig.  4  shows  the  rdf  as  a  function  of  radial  distance  r  for  all 
residues  and  phosphates  at  t  =  106  in  simulations  where  CK sl 
=  0.03  and  0.05.  For  each  eKph  in  simulations  with  phosphate, 
it  is  apparent  that  residues  are  concentrated  at  set  intervals  as 
the  radial  distance  is  increased  from  each  residue  (A).  This 
suggests  that  peptides  are  uniformly  patterned  within  the 
composite.  A  similar  distribution  is  also  observed  when  the 
rdf  of  phosphate  is  plotted  as  a  function  of  r  (B).  When  the  rdf 
of  residues  and  phosphates  for  each  r  are  compared,  the  results 
demonstrate  that  the  aggregates  are  organized  as  alternating 
phosphates  and  peptides.  For  each  concentrated  distribution 
of  residues  at  a  radial  distance  r,  there  is  a  corresponding  drop 
in  phosphate  groups  and  vice  versa.  This  alternating  trend  is 
observed  in  simulations  through  the  entire  range  of  eKph-  The 
findings  confirm  that  the  aggregates  are  formed  through 
peptide  molecules  interconnected  via  phosphate  bridges.  The 
patterning  of  peptides  within  the  aggregates  is  mediated  by 
phosphate,  because  this  unique  distribution  of  residues  is  not 
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Fig.  4  The  rdf  for  all  residues  (A)  and  for  phosphates  (B)  for 
simulations  where  the  peptide  concentration  (CKSL)  is  0.03  (left)  and 
0.05  (right).  The  legend  (inset)  lists  eKph-  The  rdf  plotted  for  all 
residues  in  simulations  without  phosphate  is  listed  in  the  legend  as  0.0. 
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observed  in  the  absence  of  phosphate  (inverted,  filled  triangles 
in  top  row  of  graphs). 

3.4  Peptide  mobility  and  conformation 

We  examined  the  mobility  of  each  residue  by  monitoring  its 
movement  throughout  the  simulation.  Residue  mobility  is 
defined  as  the  number  of  successful  moves  per  unit  time  step 
MCS.  Fig.  5  shows  the  mobility  profile  of  KSL  by  measuring 
the  number  of  successful  moves  of  each  residue  in  the  absence 
and  presence  of  phosphate  and  over  each  eKph-  The  mobility  of 
residues  is  much  higher  in  the  absence  of  phosphate.  The 
residues  at  the  end  of  the  peptide  chains  (i.e.,  Ki  and  Ki0)  are 
much  more  mobile  than  the  interior  residues  as  one  would 
expect  from  the  segmental  dynamics  of  a  chain.  The  mobility 
of  interior  residues  is  nearly  three  times  lower  than  the  end 
residues  but  still  much  higher  than  the  corresponding  residues 
in  the  presence  of  phosphate.  A  different  mobility  profile  is 
recorded  for  peptide  in  the  presence  of  phosphate.  In  general, 
the  mobility  of  each  lysine  is  lower  relative  to  the  hydrophobic 
residues  mobility.  We  see  that  the  mobility  of  each  lysine 
residue  decreases  monotonically  on  increasing  eKph-  The 
mobility  of  lysine  decreases  over  increasing  eKph  due  to  its 
selective  binding  with  the  phosphate.  Phosphate  is  thus  critical 
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Fig.  5  The  average  mobility  of  each  residue  in  the  absence  of 
phosphate  (A)  and  in  the  presence  of  phosphate  (B).  In  (A),  CKsl  in 
simulations  are  0.03  (top  graph)  and  0.05  (bottom  graph).  In  (B),  the 
CKsl  are  0.03  (left)  and  0.05  (right)  and  eKph  for  each  bar  graph  is 
shown  in  parentheses.  Each  residue  is  represented  by  a  number  (shown 
at  bottom). 


in  orchestrating  the  mobility  of  each  residue  and  therefore 
attributes  its  structural  characteristics  in  a  specific  sequence 
of  KSL. 

We  also  examined  the  conformation  of  individual  peptide 
chains  as  they  self-assemble  in  the  absence  and  presence  of 
phosphate  by  calculating  the  radius  of  gyration  ( Rg )  and 
end-to-end  distance  ( Re )  of  peptide  chains  over  time  (Fig.  6). 
It  is  apparent  that  the  peptide  chains  remain  stretched  out  in 
the  absence  of  phosphate  and  contract  when  phosphate  is 
added  to  the  simulation.  A  large  change  in  peptide  conformation 
occurs  if  the  interaction  strength  between  lysine  and  phosphate 
is  changed  from  —0.5  to  —0.6.  At  eKph  =  —0.5,  peptides  slowly 
condense  into  small  structures  and  do  not  reach  an  equilibrium 
at  the  last  time  step,  while  in  simulations  where  eKph  >  —0.6, 
peptides  in  the  presence  of  phosphate  form  an  energy- 
minimized  conformation  more  rapidly.  The  final  energy- 
minimized  conformation  is  about  one-third  smaller  than  the 
peptide  in  its  fully  extended  and  freely  mobile  state. 

A  reduction  in  Rg  and  Re  further  demonstrates  that  the 
peptides  aggregate  and  are  either  compacted  by  lysine- 
phosphate  interactions  or  fold  upon  themselves  via  multiple 
interactions  between  a  single  phosphate  and  more  than  one 
lysine  from  the  same  peptide.  Both  of  these  interactions  are 
observed  in  final  snapshots  (Fig.  1  and  Fig.  S2  of  the  ESI|).  A 
reduction  in  Rg  and  Re  did  not  occur  in  the  absence  of 
phosphate,  even  though  we  did  observe  that  a  few  peptides 
are  folded  and  condensed  in  final  snapshots  of  these  simulations. 
Apparently,  these  types  of  peptide  conformations  are  not 
stable  in  the  absence  of  phosphate,  and  repulsive  interactions 
between  lysines  and  between  hydrophilic  and  hydrophobic 
residues  compete  with  attractive  hydrophobic  interactions  to 
maintain  random  conformations  of  mostly  extended  peptides 
throughout  the  simulation. 

4.  Summary  and  conclusions 

A  coarse-grained  model  was  used  to  study  the  self-organizing 
behavior  of  an  antimicrobial  peptide  (KSL)  and  phosphate. 
Specificity  of  each  residue  and  their  interactions  with  each 
other  and  phosphate  was  considered  in  simulations.  Each 


Fig.  6  Peptide  chain  radius  of  gyration  (A)  and  end-to-end  distance 
(B)  for  CKsl  =  0-03  (top)  and  CKsl  =  0.05  (bottom).  The  interaction 
strength  of  each  simulation  is  listed  in  the  legend  (inset).  Simulations  in 
the  absence  of  phosphate  are  listed  as  0.0. 
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component  was  allowed  to  execute  its  stochastic  motion  via 
the  Metropolis  algorithm.  Both  local  and  global  physical 
quantities  were  monitored  during  the  course  of  simulations. 
The  local  quantities  were  mobility  profiles  of  each  residue  and 
their  characteristic  neighborhood  within  the  range  of  inter¬ 
action.  The  global  quantities  included  the  radial  distribution 
function  {rdf),  the  radius  of  gyration  ( Rg )  and  end-to-end 
distance  ( Re )  of  peptide  chains. 

We  find  that  peptide  assembly  into  scaffold  structures  is 
mediated  by  phosphate.  Upon  assembly,  the  peptides  and 
phosphates  segregate  and  form  composite  aggregates.  The 
simulations  reveal  that  when  assuming  stronger  interactions 
between  peptide  and  phosphate,  the  mobility  of  lysine  residues 
decreases,  while  the  average  number  of  surrounding  phosphate 
within  the  range  of  interactions  increases.  Thermodynamic 
equilibration  directs  phosphate  constituents  to  mobilize 
around  lysine  residues,  recruiting  peptides  to  form  composite 
aggregates  involving  scaffolds.  Monitoring  the  Rg  and  Re  of 
peptides  over  time  in  simulations  where  eKph  >  —0.5  reveals 
that  the  local  peptide  structure  condenses  relatively  quickly. 
At  higher  eKph,  most  of  the  phosphates  become  strongly 
mobilized  around  lysine,  resulting  in  a  compact  conformation 
of  the  peptide.  As  shown  by  the  antagonistic  fluctuations  of 
the  residue  and  phosphate  rdf  as  a  function  of  radial  distance 
r,  the  equilibrated  state  is  an  ordered  assembly  of  peptide 
residues  interspersed  with  phosphate. 

The  specific  ordered  arrangement  of  phosphate  and  KSL  is 
an  intriguing  result  and  reflects  several  proposed  mechanisms 
addressing  the  molecular  self-assembly  of  biosilicifying 
polymers.  Phosphate-bridging  mechanisms  have  been  proposed 
to  explain  the  aggregation  of  peptides  and  polyamine  chains 
prior  to  inducing  silica  formation.6,7,11,34,35  The  polyionic 
structure  of  phosphate  facilitates  bonding  to  multiple  peptides 
and  the  cross-linking  required  to  build  supramolecular  structures 
of  aggregated  peptide.  Phosphate  cross-linking  of  biomolecules 
prior  to  silica  condensation  has  been  shown  in  poly-L-lysine 
supramolecular  assemblies,34,35  poly(allylamine)  aggregates,11 
and  naturally-occurring  LCPA  scaffolds.36  In  the  bio¬ 
mineralization  of  silica,  it  is  hypothesized  that  the  coordination 
of  phosphate  and  peptide  may  induce  a  conformational  shift 
to  spatially  or  electrostatically  orient  functional  groups  to  bind 
and  react  with  silicate  species.  This  has  been  proposed  in 
silicification  reactions  using  non-peptidyl  polyamines13,37,38 
and  may  likely  occur  in  similar  reactions  with  amine-rich 
polypeptides.  The  results  of  our  simulations,  which  demon¬ 
strate  a  highly  ordered  assembly  of  peptide  in  the  presence  of 
phosphate,  are  consistent  with  proposed  models  of  supra¬ 
molecular  polyamine  and  peptide  assemblies  that  act  as  silica 
template. 

The  phosphate  bridging  interactions  observed  in  our 
simulations  also  reflect  the  proposed  mechanisms  of  in  vitro 
silaffin  self-assembly  prior  to  silica  deposition.6,7  The 
self-assembly  of  native  silaffin  peptide  fragments  is  mediated 
by  intermolecular  interactions  between  the  cationic  polyamine 
chains  on  lysines  and  anionic  phosphoryl  modifications  on 
serines.  When  the  peptide  is  synthesized  without  any  modifications 
(SSKKSGS Y SGSKGSKRRIL,  known  as  the  R5  peptide), 
assembly  into  aggregates  is  induced  in  the  presence  of  phos¬ 
phate  buffer.  Without  the  phosphoryl  groups  on  the  serines  of 


the  R5  peptide,  phosphate  buffer  is  required  for  self-assembly. 
If  the  serines  are  present  in  the  R5  peptide  sequence  solely  to 
tether  phosphate  molecules  and  the  intrinsic  properties  of  the 
serine  residues  are  not  directly  involved  in  self-assembly,  the 
remaining  sequence  of  R5  (KKGYGKGKRRIL)  shares 
significant  identity  and  arrangement  to  KSL,  except  for  the 
C-terminal  RRIL  residues.  While  the  RRIL  motif  has  shown 
to  control  assembly  of  R5  into  large  aggregates,  it  is  not 
essential,  because  the  same  peptide  without  the  RRIL  motif 
(SSKKSGS  YSGSK)  also  induces  the  formation  of  active  silica 
precipitating  assemblies.39  Considering  the  sequence  of 
cationic  and  hydrophobic  residues  on  the  R5  peptide  is  nearly 
identical  to  KSL,  it  would  not  be  surprising  if  these  two 
peptides  share  similar  assembly  mechanisms.  Similar  modeling 
studies  using  the  R5  peptide  are  underway  and  comparison 
between  R5  and  KSL  may  assist  in  a  better  understanding  of 
the  fundamental  mechanisms  that  drive  peptide  self-assembly. 

Our  coarse-grained  model  is  able  to  characterize  the 
structural  properties  of  peptide  assemblies  that  result  from 
residue-specific  interactions  and  also  global  dynamics  of  the 
system,  because  of  the  assimilation  of  the  intra-residue  details 
in  a  single  peptide  node.  In  contrast,  simulations  that  define 
the  atomic  structure  of  each  residue  of  KSL  and  phosphate 
would  limit  the  number  of  constituents  and  time  scales 
necessary  for  ordered  aggregates  to  form  and  therefore  cannot 
accurately  define  all  structural  dynamics  of  the  complete 
system.  Though  to  identify  the  molecular  structure  and 
electronic  state  of  peptide-phosphate  templates  that  induce 
silica  condensation,  atomistic  simulations  are  required  to 
characterize  the  specific  coordination  chemistry  between 
functional  groups  of  amino  acids  and  phosphate.  Subsequent 
work  will  build  upon  the  simulation  model  presented  here  and 
more  accurately  define  the  atomic  structures  of  KSL  and 
phosphate  and  also  survey  additional  biomineralizing  peptides 
in  different  ionic  solvents  to  reveal  how  different  sequences  of 
amino  acids  interact  to  form  specific  supramolecular  assemblies. 
Once  a  repertoire  of  interactions  has  been  identified,  we  can 
then  move  forward  to  design  molecular  blueprints  of  complex 
structures  that  will  mimic  the  vast  diversity  of  nanostructured 
materials  we  find  in  nature. 
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